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ABSTRACT 

We present the first and so far the only simulations to follow the fine-grained phase-space 
structure of galaxy haloes formed from generic ACDM initial conditions. We integrate the 
geodesic deviation equation in tandem with the N-body equations of motion, demonstrating 
that this can produce numerically converged results for the properties of fine-grained phase- 
space streams and their associated caustics, even in the inner regions of haloes. Our effective 
resolution for such structures is many orders of magnitude better than achieved by conven- 
tional techniques on even the largest simulations. We apply these methods to the six Milky 
Way-mass haloes of the Aquarius Project. At 8 kpc from halo centre a typical point intersects 
about 10 14 streams with a very broad range of individual densities; the ~ 10 6 most massive 
streams contribute about half of the local dark matter density. As a result, the velocity dis- 
tribution of dark matter particles should be very smooth with the most massive fine-grained 
stream contributing about 0.1% of the total signal. Dark matter particles at this radius have 
typically passed 200 caustics since the Big Bang, with a 5 to 95% range of 50 to 500. Such 
caustic counts are a measure of the total amount of dynamical mixing and are very robustly 
determined by our technique. The peak densities on present-day caustics in the inner halo 
almost all lie well below the mean local dark matter density. As a result caustics provide a 
negligible boost (< 0.1%) to the predicted local dark matter annihilation rate. The effective 
boost is larger in the outer halo but never exceeds about 10%. Thus fine-grained streams and 
their associated caustics have no effect on the detectability of dark matter, either directly in 
Earth-bound laboratories, or indirectly through annihilation radiation, with the exception that 
resonant cavity experiments searching for axions may see the most massive local fine-grained 
streams because of their extreme localisation in energy/momentum space. 
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1 INTRODUCTION 

Dark matter is supposed to be the principal driver of structure for- 
mation in the Universe, and a whole industry has developed over 
the last few decades searching for the still elusive dark matter 
particle. Particle physics has provided some well-motivated candi- 
dates. Among them, the neutralino, a weakly interacting massive 
particle (WIMP) associated with a supersymmetric extension of 
the standard model of particle physics, is currently favour ed (see 
ljungman et aL1ll996l : iBertone et alj|2003 ; lBergstr6mll2009l for re- 
views). Although this particle interacts with standard model parti- 
cles only weakly, it may nevertheless be detectable. Recently the 
dark matter community has been stimulated by a variety of ob- 
served "anomalies" in cosmic-ray signals. Among these are: (i) 
an excess in the posi tron fraction recent ly (re)measured by the 
PAMELA experiment jAdriani et alj2009h : (ii) an excess in the to- 
tal flux of electrons and positrons measured by the FERMI satellite 
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dAbdo et al.ll2009h : (iii) an excess in diffuse microwave radiation 
in the general direc tion of the Galactic Centre (the "WMAP haze", 
iHooper et ai]|2007h . Although these results may well find their ex- 
planation in or d inary astrophysical processes (see, for example, 
iMalvshev et al.l |2009) for a pulsar explanation of the PAMELA 
results), it is alluring to relate them all to dark matter annihilation 
in the Galactic halo. 

In addition to these indirect signals, the DAM A/LIBRA dark 
matter detec tion experiment has, for about a decade dBernabei et al.l 
I2008L |2010|) seen a significant amplitude modulation in its detec- 
tor count rate. Although a pparently incompatible with results from 
other experiments (see e.g.lSavage et al . 2004; Go ndolo & Gelmiml 
120051 : lGelminill2006l: iFinkbeiner et alj|2009h . the DAMA/LIBRA 
collaboration interprets this as an annual modulation of the flux of 
dark matter particles through their detectors. 

The small-scale structure of Cold Dark Matter (CDM) haloes 
can, in principle, substantially influence the signal both in indirect 
and in direct dark matter detection experiments. Annihilation rates 
and detector count rates depend strongly on the density and en- 
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ergy distributions of particles local to other particles and to detec- 
tor nuclei, respectively, and so are sensitive to small-scale fluctu- 
ations in these distributions. For example, for cross-sections con- 
sistent with relic abundance constraints, the PAMELA data can be 
explained through annihilation only if rates are 100 to 1000 times 
those predicted for a locally smooth dark matter distribution, thu s 
requiring substantial "boost factors" (e.g. iBergstrom et alj|2008l) . 
Such boosts could come from non-standard particle physics, or they 
might reflect substantial inhomogeneities in the dark matter distri- 
bution on small scale. A population of abundant, self-bound sub- 
haloes of very low mass but very high interna l density has often 
been invoked in this context l lMooreetal .11 19991 ; iGreen et al .1120051 ; 
iDiemand et al ]|2005l . l2008l) . but recent high-resolution simulations 
suggest that such subhaloes are neither dense enough nor abun- 
dant enough in the inner regions of ACDM haloes to have more 
than a minor effect o n the observable signatures of annihilation 
dSpringel et"al]|2008bh . 

Another possible boost mechanism is related to the fine- 
grained structure of dark matter haloes. Before the onset of non- 
linear structure formation, dark matter was almost uniformly dis- 
tributed, with weak density and velocity perturbations and with 
very small thermal motions. The particles thus occupied a thin, 
space-filling and almost three-dimensional sheet in the full six- 
dimensional phase-space. Subsequent collisionless evolution under 
gravity stretched and folded this sheet, but did not tear it. Thus, the 
dark matter distribution at a typical point in a present-day halo is 
predicted to be a superposition of many fine-grained streams, each 
of which has a very small velocity dispersion, has a density and 
mean velocity which vary smoothly with position, and corresponds 
to material from the vicinity of a different point in the linear ini- 
tial conditions. Folds in the fine-grained phase-sheet give rise to 
projective catastrophes known as caustics where the spatial density 
and hence the annihilation rate is locally very high, limit ed only 
by the small but nonzero thickne ss of the phase-sheet (e.g. lHoganl 
1200 ll ; iNataraian & Sikividl2008t> . Although the extraordinary im- 
provement of N-body simulations in recent years has allowed many 
aspects of the dark matter distribution at the solar position to be 
predic t ed in considerable detail (see, for exam ple, ISpringel et al.l 
l2008al : IDiemand et aTll2008l : IStadel et alj|2009h , such simulations 
are still very far from resolving the fine-grained phase-space struc- 
ture which gives rise to caustics. They are thus unable to provide 
a realistic estimate of how much caustics boost annihilation in the 
inner Galactic halo. 

Fine-grained structure might also play a crucial role in the 
interpretation and modelling of dark matter signals in laboratory 
detectors like DAMA/LIBRA. It is unclear, for example, whether 
the Maxwellian usually assumed describes the dark matter velocity 
distribution at the solar position accurately. Indeed, recent simula- 
tion work has shown that a multivariate Gaussian is likely a better 
description, and that the assembly history of the Milky Way may 
be reflected in broad feat ures in the particle energy distribution 
dVogelsberger et alj|2009"3) . One may also wonder whether indi- 
vidual fine-grained streams might eventually be visible in detector 
signals as spikes at specific velocities. Axion detectors like ADMX, 
in particular, have very high energy resolution, and could in princi- 
ple disentangle hundreds of thousands of streams. Again this is far 
beyond the resolution limit of current analysis techniques applied 
to even the highest resolution N-body simulations of halo forma- 
tion. 

In recent work we have developed an entirely new approach 
capable of following the evolution of fine-grained streams and 
their associated caustics in fully general simulations of halo for- 



mation. This approach is based on integrating the geodesic de- 
viation equation (GDE) in tandem with the N-body equations of 
motion. For every simulation particle and at every timestep it 
returns the spatial density and the velocity dispersion tensor of 
the fine-grained stream in which the particle is embedded. This 
in turn allows all passages o f the particle through a c austic to 
be identified and recorded. In [ Vogelsberger et al.l | |200§) we pre- 
sented this GDE scheme and tested it on orbits in fixed poten- 
tials and on N-body simu lations of static equilibrium haloes. In 
I Vogelsberger et alj j2009bh we then applied the scheme to a simpli- 
fied model of CDM halo formation: collapse from spherically sym- 
metric, self-similar and linear initial conditions. This showed that 
earli er work based on spherica l ly symmetric similarit y solutions 
(e.g. INataraian &Sikiviell2006l : iMohavaee et alll2007t) was unre- 
alistic because these solutions are violently unstable to nonradial 
perturbations. These substantially alter t he stre am and caustic struc- 
ture. Finally in lWhite & Vogelsbergerl ( 120091) we showed how the 
GDE scheme allows the annihilation enhancement due to caustics 
and other fine-grained structure to be calculated explicitly by time- 
integration along particle orbits rather than by spatial integration 
over the particle distribution. The current paper is the culmination 
of this programme and analyses for the first time the fine-grained 
structure of haloes forming from fully general ACDM initial condi- 
tions. We resimulate haloes from the Aquarius Project, which stud- 
ied six Milky Way mass objects at various numerical resolutions. 
The coarse-grained structure of these haloes has already been stud- 
ied in con siderable detail in prev ious papers, e.g. the s ubhalo pop- 
ulation in ISpringel etall d2008al) and ISpringel et all d2008bh, the 
dark m atter distribution in the inner regions in I Vogelsberger et all 
d2009al) and various radial profiles in lNavarro et alJd2010l) . 

The plan of our paper is as follows. In Section 2 we briefly 
describe the initial conditions and the numerical approach we use 
to resolve the fine-grained phase-space structure. Section 3 begins 
by presenting results on one of the Aquarius haloes at a variety of 
resolutions. We demonstrate that, for fixed gravitational softening, 
convergent results can be obtained, free of significant discreteness 
or two-body relaxation effects. We then analyse the implications 
for dark matter detection of the structure predicted for fine-grained 
streams and caustics. Finally, we use our full halo sample to tackle 
the question of how much scatter in fine-grained properties is ex- 
pected. We make some concluding remarks in Section 4. An Ap- 
pendix shows that although discreteness effects are negligible in 
our high resolution simulations, their predictions for fine-grained 
stream densities (but not for caustic structure) are influenced by 
gravitational softening because this affects the tidal forces on par- 
ticles passing close to halo centre. 



2 INITIAL CONDITIONS AND NUMERICAL METHODS 

We resi mulate the six Milky Way-mass haloes of the Aquarius 
Project dSpringel et all l2008al) using our GDE technique to fol- 
low the fine-grained phase-space evolution in detail. The cos- 
mological parameters assumed for these ACDM simulations are 
Q m0 = 0.25, Qao = 0.75, as = 0.9, n s = 1 and H = 
73 km s" 1 Mpc -1 , where all quantities have their standard def- 
initions. The haloes were selected at z — based on their mass, 
and were required to have no close massive companion at that time; 
they are named Aq-A to Aq-F. Each halo was resimulated at a va- 
riety of numerical resolutions, indicated by a number n in its full 
name, e.g. Aq-A-3. The resolution levels differ in the particle mass 
and softening length employed, with 1 designating the highest res- 
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olution and 5 the lowest. In the following we will use the same 
naming convention and mass resolution as in the original papers 
but different softening lengths. This is because accurate integration 
of the geodesic deviation equation (GDE) requires a larger soft- 
ening length to achieve stable results than does integration of the 
particle trajectories themselves. Unless otherwise stated, we use a 
constant comoving Plummer-equivalent softening of 3.4 kpc in all 
our simulations. The Appendix examines the effect of softening on 
our results in some detail . All ot her simulation parameters are the 
same as in lSpringel et alj d2008ai) . to which we refer the reader for 
further technical information. 

Our experiment s are carried out with the P-GADGET-3 
code dSpringell 1200^ with the GDE modifications described in 
IVogelsberger et alj J2008h and IVogelsberger et alj d2009bh . where 
we focused on a description of the relevant equations in physical 
coordinates and applied them to static potentials, to isolated equi- 
librium haloes and to halo formation from self-similar, spherically 
symmetric initial conditions. The simulations of this paper take into 
account the full ACDM framework, using the cosmological param- 
eters given above, and are carried out in comoving coordinates. We 
therefore begin by describing how the GDE and the associated ten- 
sors can be transformed from physical to comoving coordinates, 
and how this is implemented in our simulation code. For complete- 
ness, we also describe how the physical stream density is calculated 
and how caustic passages can be identified. 

The time-dependent transformation to the comoving frame 
used by our simulation code is given bjQ 
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Figure 1. Spherically averaged density profiles of Aq-A-3, 4 and 5 at z = 0. 
The comoving softening length of all simulations is the same (e = 3.4 kpc; 
note that T200 = 246 kpc for this halo) but the particle number increases 
from Aq-A-5 to Aq-A-3 as shown by the values of A^oo in the figure. The 
lo wer panel shows the d ifference between these profiles and that presented 
bv lNavarro et alj 1201 (J) for Aq-A-1, for which Af 2 oo ~ 10 9 . The conver- 
gence between the different resolutions is very good. All begin to fall below 
Aq-A-1 for r < 0.025r2oo = 6 kpc because of the larger softening of our 
GDE resimulations. An arrow in the lower panel indicates 8 kpc, the radial 
position of the Sun within the Milky Way. 



where H — a/a denotes the Hubble parameter, a the scale factor, 
and comoving coordinates are primed. In physical coordinates the 
distortion tensor is defined as 



dq 



dx/dq dx/dp 
dv/dq dv/dp 
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where q and p denote the initial position and velocity (in physical 
coordinates), respectively. It describes how a local phase-space el- 
ement in physical coordinates is deformed along the trajectory of 
the particle while conserving its volume in order to obey Liouville's 
theorem. In comoving coordinates the distortion tensor is accord- 
ingly defined as 
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where all phase-space coordinates are now expressed in comoving 
coordinates. The relation between physical and comoving distor- 
tion tensors can be derived from the differential relations, yielding 



D(t) = CV^(t)L> (t)CW (tinitial), 

where we have defined the two transformation tensors 
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1 We use the following notation to clearly distinguish between three- and 
six-dimensional quantities: an underline denotes a R 3 vector and two of 
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of them denote a R 6x6 m atrix. See also Vogelsberg er et al] J2008h and 
Vogelsberg er et alj l2009bl) . 
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We note that one of these transformation tensors is evaluated at 
t whereas the other is evaluated at the initial time imitiai- The 
reason for this is the time-dependence of the coordinate trans- 
formation in Eq. 1, where both the scale factor and the Hubble 
parameter change with time. Liouville's theorem guarantees that 

det(D(t)) = det(D (t)) — 1, so the comoving distortion ten- 
sor has the same conserved determinant as the physical one. This 
means that local phase-space elements in the comoving frame con- 
serve both volume and orientation as the system evolves. 

To integrate the comoving distortion tensor, we need to derive 
its equations of motion. We start with the physical GDE equation 

15 = TTj, (7) 
and introduce the peculiar potential field 

= a$ + ^r' 2 , (8) 
which is related to the density field via Poisson's equation 

Yl t </> = 47rG(p'(w')-p' b ), (9) 

where p' b denotes the comoving mean background density and the 
Laplacian is taken in comoving space. We note that the force field 
driving the motion of simulation particles is also derived from this 
peculiar potential. Since the GDE is directly related to the equa- 
tions of motion of the particles themselves, it is natural to intro- 
duce a peculiar configuration-space tidal tensor T' with compo- 
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nents TL = —d 2 cj>/dx' i dx'j 



\t' = t-^i. 



(10) 



We can use this tidal tensor to write the equations of motion for the 
comoving distortion tensor in the following form 



d'., = Xd' „ b',=±D', 



(11) 



D = -T'D' ,, D , 

=v'q' a ==x'q' =v'p' a ==x'p' 



-T'D' 



Thus, the equations of motion for the comoving distortion tensor 
have exactly the same form as those for the physical one. The 
only difference is the appearance of the scale factor a in the co- 
moving equations. The initial conditions for these equations follow 
from those for the physical distortion tensor, taking into account 
T'-i - 71- -, — T 

^ X — >X ^ X — VX ^ 



If , , (tinitial) = i, ry , , (t in itial) = 0, 

^=x q v = ^=x'p' K ' — 

D' , , (tinitial) =0, t>' , , (tinitial) = 1. 



(12) 



Based on these equations we can construct kick- and drift-operators 
for the leapfrog time integrator, similar to those for the position and 
velocity in the comoving frame. 

Let us now discuss how to derive stream densities from these 
comoving quantities. Each stream density is associated with a 
particular simulation particle, and describes the local density of 
the fine-grained stream it is embedded in. Our goal is to derive 
an equation that yields the physical stream density from the co- 
moving distortion tensor discussed above. To begin, we note that 
the configuration-space distortion tensor can be derived from the 
phase-space distortion tensor by applying two projection operators 
as follows 



1)D\ = 

=<i 



(13) 



where V = dV_/dq is the spatial gradient of V_(q), the 
mean initial sheet velocity as a function of initial position (see 
IWhite & Vogelsbergedl2009l. for a more formal definition). Caus- 
tic passages can be identified by sign changes in the determinant 
of this tensor. We note that this is true both in physical and in co- 
moving space. The physical stream density can be calculated from 
the volume stretch factor implied by this linear transformation in 
physical configuration-space 



Ps,0 
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as we described in IVogelsberger et all d2008i ). Using p s .o = 
Ps,o/°(£initiai) 3 and defining the peculiar velocity gradient V_' t = 
V_ — //(tinitial)! we can write the physical stream density as: 



Ps,0 



\det\D', , + a(t initia i) 2 V! , r/ , II 

|_ x' q' q' x'p'J j 



(15) 



At early times when the dark matter is nearly uniform, all velocities 
and thus all velocity gradients are small and the second term in the 
determinant in this equation is small compared to the first. Thus we 
may approximate 



Ps = 



Ps,0 



(lot 



(16) 



This equation is sufficiently accurate for our purposes and we 
will use it below to calculate the stream densities associated with 
each particle in our simulations. For consistency we will also ne- 
glect the (small) initial variations in density and will assume the 
cosmological mean value everywhere, i.e. 



Ps,0 



8ttG 



(17) 



3 QUANTIFYING FINE-GRAINED HALO STRUCTURE 

The results in the next few subsections are based on resimulations 
of the Aq-A halo at resolution levels 3, 4 and 5. We note again 
that our simulations use a significantly larger softening length than 
the original Aquarius simulations in order to mitigate the numeri- 
cal noise sensitivity of the GDE. Unless otherwise stated, we use 
a Plummer-equivalent comoving softening length of 3.4 kpc. Note 
that we do not change softening between the different resolution 
levels that we are going to discuss. We will demonstrate that we 
reach convergence at level 4 with this set-up. Additional experi- 
ments show that this is the smallest softening length for which we 
could achieve convergence at this resolution level; smaller values 
require a larger particle number to converge. In the Appendix we 
show that our results for fine-grained stream densities (but not for 
caustic counts) are sensitive to the assumed softening because of 
its effect on the central structure of halo and subhalo potentials. 
Where needed, we also assume a dark matter particle with a veloc- 
ity dispersion of 0.03 cm/s today. This corresponds to a neutralino 
of mass 100 GeV/c 2 that decoupled kinetically at a temperature 
around 10 MeV. 

To get started we show in Fig.Q]the spherically averaged den- 
sity profile of the Aq-A halo at redshift 2 = for all three reso- 
lution levels. Note that 7-200 = 246 kpc for this halo. Convergence 
is excellent over the full radial range plotted, as is more obvious in 
the lower panel w here we plot the diffe rence between our profiles 
and that given by iNavarro et al. I j20ld) for the highest resolution 
simulation Aq-A- 1 . Softening clearly effects all three of our simu- 
lations similarly and only in the innermost regions; the deviation is 
at the percent level at the radius corresponding to the Sun's position 
within the Milky Way, but reaches 10% by 4 kpc. INavarro et al .1 
d2010l) made similar convergence tests for the radial profiles of a 
variety of quantities in the original simulations. 



3.1 Caustic counts 

Our focus in this paper is not on the coarse-grained but on the fine- 
grained structure of our haloes. We begin by looking at the dis- 
tribution of the number of caustics passed by particles in Aq-A. 
Caustics are identified through changes in the sign of the determi- 
nant of the comoving configuration-space distortion tensor and are 
counted along the trajectory of each particle. In Fig. [2] (top panel) 
we show the distribution of this caustic count in a phase-space dia- 
gram at z — for Aq-A-3. Colour encodes the number of caustics 
passed by each particle, as indicated by the colour bar. The parti- 
cles were sorted by their caustic count before plotting, so that the 
particle with the highest caustic count is plotted in each occupied 
pixel. This phase-space diagram c an be compared directly to the 
one in IVogelsberger et alj J2009bh which shows an isolated halo 
which grew from a spherically symmetric and radially self-similar 
perturbation of an Einstein-De Sitter universe. In both cases the 
number of caustic passages increases towards the centre of the halo, 
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0.5 1.0 1.5 2.0 2.5 

r / r 200 

Figure 2. Top panel: Phase-space structure of the Aq-A-3 halo at z = 0. Colour encodes the number of caustics passed by each individual simulation particle 
with the largest numbers plotted last. Subhaloes clearly stand out in this plot, because shorter orbital times lead to a higher number of caustic passages in 
their centres. We note that the streams visible in this plot are not fine-grained streams, but rather tidal streams resulting from disrupted subhaloes. Particles 
previously associated with subhaloes have higher caustic counts than other neighboring main halo particles, so they stand out in this phase-space portrait. 
Bottom panel: Same as the top panel, but with all bound subhaloes removed so that only the main halo component remains. Tidal streams from disrupted 
subhaloes are now more clearly visible. Note that this plot only shows particles bound to the main halo, while the top panel includes all particles within 5 T200 ■ 
This is why the main halo contribution ends at about 2.5 r200- The colour scales for the two panels are the same, as indicated in the top panel. 
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Figure 4. Slices of thickness 0.1 Mpc through (xy, xz, yz) projections of the particle distribution in Aq-A-5, filtered by the number of caustics passed. The 
first four rows match the exact caustic passage number indicated, whereas the last row filters for particles that have passed more than 50 caustics. The degree 
of structure in the different panels increases with the number of caustics passed. Particles that have not passed any caustic form a smooth distribution compared 
to the other panels of the figure. 
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Figure 5. Median number of caustic passages (solid thick lines) in Aq-A 
as a function of radius. In the outer regions subhaloes show up as peaks. 
The convergence between the different resolutions is very good. Thin lines 
show the upper and lower 25%, 5% and 1% quantiles of the caustic count 
distribution in each radial bin for Aq-A-3 and Aq-A-4. The agreement be- 
tween the distributions is excellent, at all radii except for some shifts at 
high counts in the outer regions which reflect small offsets in the positions 
of substructures. Clearly, caustic identification is stable against discreteness 
effects at these resolutions. 



due to the shorter dynamical timescales in the inner regions: caus- 
tic count is roughly proportional to the number of orbits executed 
by each particle, because caustics occur near orbital turning points. 
The most important difference between the isolated and ACDM 
haloes is in the subhaloes and associated tidal streams seen in the 
ACDM case. Subhaloes stand out clearly in Fig.[2]since their parti- 
cles have short orbital periods and undergo many caustic passages. 

In the bottom panel of Fig. [2] we eliminate self-bound sub- 
haloes to show only particles belonging to the main halo. Here tidal 
streams corresponding to disrupted subhaloes can be traced very 
well. Material from such disrupted subhaloes has a different caustic 
count distribution than other main halo material, and so stands out 
in these phase-space plots. In Fig. [3] we show two zooms into the 
phase-space structure around the most prominent subhalo of Fig. [2] 
These show the structure of the tidal streams in more detail. We 
note that these tidal streams are not equivalent to the fine-grained 
streams we will discuss below. As we will see, the latter are not 
visible in such phase-space plots because of their large number and 
their low densities. 

Caustic counts can be used to study how dynamical mixing is 
related to particle location. To this end, we filter the particle dis- 
tribution by number of caustic passages and plot the result in three 
orthogonal projections in Fig. [4] The top row of this figure shows 
that particles that have passed no caustic are almost homogeneously 
distributed, with no clear structure. Particles that have passed ex- 
actly one caustic already delineate the large-scale structure of the 
density field. The main filament passing through the primary halo is 
visible, as is the halo itself and additional smaller filaments point- 
ing towards it. The halo is even more dominant in the distribution 
of particles that have passed two or ten caustics. The bottom row of 
Fig. E] shows particles that have passed at least 50 caustics. These 
are found only near the centres of the main halo and of the more 
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Figure 6. Top panel: Distribution of number of caustic passages at z = 
for particles in a set of nested spherical shells in Aq-A-3 for our standard 
softening of 3.4 kpc. As expected, higher caustic counts are found in the 
inner regions. Bottom panel: Caustic count distributions within 160 kpc for 
Aq-A-3 resimulations with various softenings. The smaller the softening 
length, the longer the tail of high-count particles. The particles affected are 
almost all in the innermost region of the the main halo. Note that aside from 
this tail, the histograms are almost independent of softening. 



massive subhaloes. These are the densest regions with the shortest 
orbital times. Typical values of the caustic count in various regions 
thus indicate their level of dynamically mixing. Large values corre- 
spond to well-mixed regions whereas small numbers correspond to 
dynamically "young" regions. Particles that have passed no caustic 
are still in the quasilinear phase of structure growth. 

We can compress the caustic count information in a radial pro- 
file. This is shown in Fig. [5] Profiles of the median caustic count 
for Aq-A-3,4 and 5 are plotted as solid curves. As already demon- 
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Figure 7. The distribution of stream density in units of the cosmic mean 
density (pi,) as a function of radius within halo Aq-A. Continuous coloured 
curves give the median stream density of particles as a function of their 
distance from halo centre for Aq-A-3 (blue), Aq-A-4 (green) and Aq-A-5 
(red). Dashed and dotted curves give the 0.5, 2.5, 10, 90, 97.5 and 99.5% 
points of the distribution of stream density at each radius for Aq-A-3 (blue) 
and Aq-A-4 (green). The stream density distributions of the two higher res- 
olution simulations agree well, but Aq-A-5 gives somewhat lower stream 
densities (for all percentiles, not just the median shown) as a result of its 
higher discreteness noise, We conclude that the resolution of Aq-A-4 is suf- 
ficient to suppress 2-body relaxation and other discreteness effects to an 
acceptably small level. Note that, remarkably, the median stream density 
is within an order of magnitude of the cosmic mean at all radii. For com- 
parison, the black solid curve repeats the mean mass density profile of the 
halo from Fig.[T] Note that although the distribution of stream density spans 
twelve orders of magnitude in the inner halo, fewer than 1 % of dark matter 
particles are in streams with densities exceeding 1 % of the local mean. 



strated in lVogelsberger et"all d2008h . the numerical identification of 
caustics is very robust, and as a result the number of caustic pas- 
sages is little affected by numerical noise. This is why the median 
profile is almost independent of resolution in Fig. [5] with remark- 
ably good agreement in the inner halo. In the outer regions sub- 
haloes show up as peaks in the caustic count and small shifts in 
their position between the different simulations show up as appar- 
ent noise. For the highest resolution simulation, Aq-A-3, we also 
plot the upper and lower 1, 5 and 25% points of the count distri- 
bution at each radius. These parallel the distribution of the median 
count and span an order of magnitude at each radius. The median 
profile can be compared to that given bv lVogelsberger et aljj2009bh 
for an isolated halo growing from self-similar initial conditions. At 
a given fraction of the virial radius, the typical number of caustic 
passages is only a few times larger in the more complex ACDM 
case. 

The increase in caustic count towards halo centre can also be 
seen in Fig. [6] (top panel), where we present count histograms for 
particles in a set of nested spherical shells. The shift of the distribu- 
tions towards lower counts with increasing radius is very obvious, 
and within about 30 kpc this is accompanied by a suppression of 
the high-count tail. At larger radii, however, the distributions al- 
ways extends out to about 200 counts. This tail is contributed by 
particles from present or recently disrupted subhaloes (see below). 
In the outer halo, the majority of particles have nevertheless passed 
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Figure 8. The number of streams in Aq-A as a function of radius, as esti- 
mated from the fine-grained stream densities of individual simulation par- 
ticles. As in Fig. [7] the Aq-A-5 results disagree with those from the two 
higher resolution simulations, but Aq-A-3 and Aq-A-4 agree reasonably 
well. Stream numbers are estimated by dividing the mean density at each ra- 
dius by an estimate of the characteristic local density of individual streams. 
For this estimate we use both the harmonic mean (solid lines) and the me- 
dian (dashed lines) stream density for particles in each radial shell. The 
harmonic mean is very sensitive to the low density tail of the stream density 
distribution, but the agreement between Aq-A-3 and Aq-A-4 remains good 
except at the smallest radii. 



relatively few caustics. As noted above, simulations using the GDE 
technique require significantly more gravitational softening than 
standard N-body simulations in order to limit discreteness noise 
in the tidal field which otherwise leads to an unphysical, nearly 
exponential decay in stream density. Caustic identification is, how- 
ever, relatively stable against such effects, so when studying caustic 
counts it is possible to reduce the softening. We show the effects of 
this in the bottom panel of Fig.[6]which compares the caustic count 
distribution within 0.16 Mpc in our standard resimulation with that 
in two additional resimulations with two and four times smaller 
softening. Smaller softening results in better resolution of the in- 
nermost regions both of the main halo and of subhaloes, and thus 
to shorter dynamical times and larger caustic counts for particles or- 
biting in these regions. This is evident as an extension of the high- 
count tail with decreasing softening. Notice, however that the re- 
gion affected is more than two orders of magnitude below the peak 
of the distribution. The shift towards higher counts only affects a 
few percent of particles that pass close to halo centre, and the dis- 
tribution is essentially unaffected below a count of about 400. 

3.2 The distribution of fine-grained stream density 

In the standard CDM paradigm, nonlinear evolution in the dark 
matter distribution can be viewed as the continual stretching, 
folding and wrapping of an almost 3-dimensional "phase-sheet" 
which initially fills configuration space almost uniformly and 
is confined near the origin in velocity space (see, for example 
IWhite & Vogelsbergej|2009h . Caustics are one generic prediction 
of such evolution. Another is that the phase-space structure near 
a typical point within a dark matter halo should consist of a su- 
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perposition of streams, each of which has extremely small veloc- 
ity dispersion and a spatial density and mean velocity which vary 
smoothly with position. If the number of streams at some point, for 
example at the position of the Sun, is relatively low, the resulting 
discreteness in the velocity distribution could give rise to measur- 
able effects in an energy-sensitive detector of the kind used in many 
dark matter experiments. To assess this possibility, we need to esti- 
mate the number and density distributions of streams at each radius. 
Our GDE formalism makes this possible by providing a value for 
the density of the fine-grained stream associated with each simula- 
tion particle. The set of stream densities corresponding to particles 
within some spherical shell is thus a (mass-weighted) Monte Carlo 
sampling of the stream densities at that radius. 

In Fig. [7] we show the distribution of fine-grained stream den- 
sity as a function of radius in halo Aq-A. Continuous curves give 
radial profiles at z = for the median stream density of particles 
at three numerical resolutions, while broken curves compare the 
tails of the distributions for the two highest resolutions. Agreement 
is excellent at all radii and throughout the distribution for Aq-A- 

3 and Aq-A-4, but the lower resolution simulation Aq-A-5 gives 
lower stream densities at all radii (and at all points in the distribu- 
tion, although we do not show this in order to avoid confusing the 
plot). This is a consequence of discreteness noise in the tidal tensor 
which affects integration of the GDE. We note that Fig.[5]and Fig.[7j 
demonstrate that the caustic count and stream density distributions 
are independent of particle number at our two highest resolutions. 
Thus, discreteness effects, in particular 2-body relaxation, have no 
significant influence on the evolution. In the Appendix we show 
that while the caustic count distribution is, in addition, indepen- 
dent of the assumed gravitational softening, the same is not true for 
the stream density distribution. This is because softening influences 
the tidal forces on particles which pass within one or two softening 
lengths of the centre of a halo or subhalo, and so influences the evo- 
lution of their stream density. As we discuss in the Appendix, this 
affects a surprisingly large fraction of halo particles at some point 
during their evolution, skewing the stream density distribution out 
to quite large radii. For our standard softening length of 3.4 kpc, 
the effects are similar in amplitude to those expected from our ne- 
glect of the baryonic component, so we do not attempt to correct 
for them. 

A remarkable result from Fig.[7jis that the median stream den- 
sity depends very little on radius, and is within an order of magni- 
tude of t he cosmic mean density at a ll radii. We found a very similar 
result in Ivbgelsberger etall d2009bl) for collapse from self-similar 
and spherically symmetric initial conditions, so it is likely to apply 
quite generally to collisionless, nonlinear collapse from a smooth 
and near-uniform early state. It implies that the dark matter dis- 
tribution from the immediate neighborhood of a randomly chosen 
point in the early universe is almost equally likely to be compressed 
or diluted relative to the cosmic mean by subsequent nonlinear evo- 
lution, and furthermore that whether it is compressed or diluted is 
almost independent of its final distance from halo centre. 

The thin dashed and dotted lines in Fig. [7] show the 0.5, 2.5, 
10, 90, 97.5 and 99.5% points of the stream-density distribution as 
a function of distance from the centre of our Aq-A-3 and Aq-A- 

4 simulations. This distribution is very broad, spanning almost 12 
orders of magnitude near halo centre. Despite this, within 0.5r2oo 
even the upper 0.5% tail of stream densities lies below the mean 
density of the halo which is shown as a solid black curve for com- 
parison. At 8 kpc, the equivalent of the Solar radius, fewer than 1% 
of all dark matter particles are part of a fine-grained stream with 
density exceeding 1% of the local mean halo density. This is a first 



10' 



7kpc<r<13kpc 
1 5kpc<r<25kpc 

35kpc<r<45kpc 

140kpc<r<180kpc 



10' 



10' 



10 : 




iog(p s /Pb. 



Figure 9. Histograms of fine-grained stream density at z = for Aq-A-3 
particles in set of spherical shells. Labels with the colour of each histogram 
give the radial range of the corresponding shell, while vertical dotted lines 
correspond to the mean halo density within that shell. All histograms are 
normalised to integrate to unity. 



Q. 

A 




10" 6 10 -4 

Ps/<P> 



Figure 10. The fraction of particles in Aq-A-3 at z = with halocentric 
radii in the range 7 to 13 kpc which have fine-grained stream density ex- 
ceeding p s is plotted as a function of p„/(p), where (p) is the mean halo 
density at 10 kpc. Half of all particles are in streams with p a > 10" 7 (p). 
See the text for a discussion of the other characteristic points marked. 
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Figure 11. The panels show a series of snapshots from the evolution of Aq-A-5. The redshift is indicated above each panel. Red and blue points refer to the 
same particles in each panel and are plotted in random order so that the probability of a random pixel being red or blue in a saturated region is proportional to 
the number of particles of each colour in that pixel. The particles were selected by separating the full z = particle distribution into a series of thin spherical 
shells centred on the main halo potential centre, and then choosing in each shell the 1% tails with the highest (blue) and lowest (red) stream density. Blue and 
red particles are thus equal in number and at z = each population is distributed in distance from halo centre in the same way as the particle distribution 
as a whole. Only particles within 10 T2oo at z = were considered. The comoving cubic region plotted lies fully within this particle set at all times and is 
centred on the centre of mass of the 200 particles which were most bound at z = 0. Clearly, low z = stream density particles typically belong to collapsed 
structures at early times, whereas high stream density particles were generally part of no structure before they were accreted smoothly onto the main halo at 
relatively late times. 



© 2010 RAS, MNRAS OOO.QJgo] 



12 M. Vogelsberger and S.D.M. White 



18.824 9.860 7.225 3.933 




x [Mpc] 

Figure 12. This figure was constructed in exactly the same way as Fig. Inl and shows the same number of red and blue particles at the same redshifts and 
within the same spatial regions. The difference is that particles were here selected on the basis of the number of caustics they have passed by z = 0, rather 
than by their present-day stream density. Thus, for the same set of thin spherical shells as before, red particles are the 1 % within each 2 = shell with the 
highest caustic count, and the blue particles are the 1 % with the lowest caustic count (random sampled if necessary among particles with equal count to obtain 
exactly the correct number). Thus the number of red and blue points in the lower right panel of this figure and of Fig. llll is the same, and they have exactly the 
same distribution in halocentric radius. The small apparent number of red points in most panels of this figure is due to the fact that high 2 = caustic count 
particles are strongly concentrated to the centres of collapsed structures at all times. The relatively small number of blue points at the earliest redshifts in this 
figure reflects the fact that many blue particles are outside the region plotted at early times. The qualitative behaviour of the distributions in this figure and in 
Fig.[TT]is similar, but the separation into high spatial density and near-uniform distributions is much more marked when particles are selected by caustic count 
rather than by stream density. In addition, low caustic count particles are accreted onto the final halo in a shell-like pattern which is a clear reflection of the 
structure seen in spherical infall models of halo formation. 
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indication that fine-grained streams are unlikely to influence direct 
detection experiments strongly. 

The total number of streams at a typical point in a radial shell 
can be estimated from the stream density distribution as the ratio 
of the mean halo density in the shell to the harmonic mean of the 
stream densities of the particles it contains. This number is dom- 
inated by the extended low-density tail of the stream-density dis- 
tribution; a very large number of very low-density streams is pre- 
dicted in the inner regions. A more relevant measure for dark mat- 
ter detection can be obtained from the median stream density of the 
particles. Half of all events in a dark matter detector will come from 
streams with densities exceeding this value, and so will come from 
a number of streams somewhat less than half the ratio of the mean 
halo density to this median value. 

In Fig.[8]we show the results of such calculations for our res- 
imulations of Aq-A at various resolutions. Solid lines show the to- 
tal number of streams at each radius estimated from the harmonic 
mean stream density, whereas dashed lines show the number of 
"massive" streams estimated from the median individual stream 
density. The two estimates differ dramatically, particularly in the 
inner halo, reflecting the very broad distribution of stream densities 
at each radius, and in particular the presence of simulation particles 
with very low stream densities. As in FigUJ convergence is some- 
what poorer than for the caustic count statistics we looked at above, 
but the results for Aq-A-3 and Aq-A-4 nevertheless agree remark- 
ably well. This is particularly notable for the solid lines, given the 
sensitivity of the harmonic mean to the low-density tail of the distri- 
bution. From Fig.[8]we estimate the total number of streams at halo- 
centric radii near that of the Sun to be around 10 14 and the number 
of "massive" streams to be about 10 6 . Clearly CDM haloes are pre- 
dicted to be very well mixed in their inner regions, and the velocity 
distribution near the Sun should appear very smooth. The very large 
number of streams predic ted near the Sun's position can be under- 
stood from the analysis in I Vogelsberger et all d2008h . which shows 
that stream density should fall inversely as the number of caus- 
tic passages in a one-dimensional system, but as the cube of this 
number for regular orbits in a three-dimensional system, and even 
faster for chaotic orbits. Given that typical particles near the Sun 
have passed a few hundred caustics, typical stream densities are 
predicted to be ~ lCP 7 /0(,. The local mass overdensity at the so- 
lar position is ~ 10 5 , so one predicts ~ 10 12 fine-grained streams 
passing through our vicinity. Additional scatter is introduced by the 
fact that a significant fraction of the particles are "pre-processed" 
in smaller mass systems before they fall into the main halo. 

We show the distribution of stream densities in a different 
form in Fig. [9] For a series of spherical shells with mean radii in- 
creasing by factors of two, we have made histograms of the fine- 
grained stream densities of Aq-A-3 particles, normalising them to 
unity so that their shapes can be compared. Beyond 30 kpc these 
histograms are all quite similar, and resemble slightly skewed log- 
normal distributions. At the bottom of our plot, three orders of mag- 
nitude below peak, they span 14 orders of magnitude in stream den- 
sity. At smaller radii, the low-stream-density tail becomes more ex- 
tended, the peak of the distribution shifts very little, and the shape 
of the high-stream-density tail is unchanged. The latter is deter- 
mined by the behaviour near caustics. As we will see below the 
maximum densities at caustics are predicted to be in the range ex- 
plored by this high-mass tail. It is notable that within 20 kpc these 
tails do not extend up to the local mean density of the halo. 

The implications of these distributions for direct detection ex- 
periments on Earth are most easily drawn from Fig. [TO] which 
shows the cumulative distribution of fine-grained stream density for 



Aq-A-3 particles with 7 kpc < r < 13 kpc at z — 0. Specifically, 
we plot the fraction of particles with stream density exceeding p s 
against p a j (p), where (p) is the mean halo density at 10 kpc. The 
fraction of particles with p a > (p) is about 2 x 10~ J , so the proba- 
bility that a single stream dominates the signal in a direct detection 
experiment (i.e. that the Earth lies sufficiently close to a sufficiently 
strong caustic) is about 2 x 10~ 5 . The fraction of particles with 
p a > 0.1 (p) is about 2 x 10~ 4 , hence the probability to see a sin- 
gle stream containing 10% of the signal is about 0.002. Similarly, 
the fraction of particles with p s > 0.01{p) is about 2 x 10 -3 , so 
a single stream containing 1 % of the signal will be seen with prob- 
ability 20%. Finally, the fraction of particles with p s > 0.001 (p) 
is about 10 -2 ; this means that at a typical "Earth" location there 
will be a few streams which individually contribute more than 0.1% 
of the local dark matter density. Thus an experiment which regis- 
ters 1000 true dark matter "events" should get a few duplicates, i.e. 
events coming from the same fine-grained stream and so having 
(very nearly) the same velocity. If the dark matter is made of ax- 
ions then a resonant detector should find an energy spectrum where 
a few tenths of a percent of the total energy density is concentrated 
in a few very narrow "spectral lines". 

3.3 Origin of the tails of the stream density and caustic count 
distributions 

Given the very broad range of stream densities and caustic counts 
present at each point within our haloes, it is interesting to inves- 
tigate how these properties are related to the dynamical history of 
individual particles. We begin by studying stream density. 

As we saw above, the "typical" stream density of particles 
(specifically, the local median value) is almost independent of local 
mean halo density, thus of local dynamical time. Lower stream den- 
sity corresponds to greater stretching of the initial phase-sheet and 
so, one might think, to more effective dynamical mixing. An an- 
ticorrelation of fine-grained phase density with local orbital time 
(i.e. with radius) is, however, at most marginally evident in the 
stream density distributions of Fig. [7] and then only in their low- 
stream-density tails. To a first approximation there is no correlation 
between stream density distribution and radius. The overall mean 
density profile of the halo is built up by having more streams of 
every density at smaller radii. Thus the relation of stream density 
to dynamical mixing is far from evident. 

We explore this issue in Fig. QT| Red and blue points here 
highlight the positions at a series of redshifts of particles which at 
z = lie within 10 r2oo of halo centre. We divided this region 
into a set of thin nested spherical shells and then for each shell we 
identified particles in the upper and lower 1% tails of the stream 
density distribution. The particles with the highest z = stream 
densities (which are typically close to caustics) are plotted blue in 
each panel while those with the lowest stream densities (which are 
typically in collapsed regions and far from caustics) are plotted red. 
Thus at z = there are the same number of red and blue points in 
the plot, and both populations are distributed in halocentric radius 
exactly as the mass distribution as a whole. The spatial distributions 
of the two populations differ markedly, however, both at z = and 
at all earlier times. Most of the red particles are part of clumps al- 
ready at early times, while most of the blue particles either appear 
diffuse or are distributed relatively smoothly through the main halo 
at all times. There is thus indeed a correlation between stream den- 
sity and dynamical history, despite the lack of correlation between 
stream density and halocentric distance in Fig. [7] 

The number of caustic passages is a clearer measure of dy- 
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Figure 13. Top panel: Radial profile of the median value of peak caustic 
density for simulations at three different resolutions. The resolution of Aq- 
A-5 is too low to get converged results, but the results from Aq-A-4 and 
Aq-A-3 agree quite well. Bottom panel: The mean halo density profile of 
Aq-A-3 is compared to the median and the quartiles of the peak density of 
caustics as a function of radius. Beyond the virial radius, typical caustics 
have peak densities which exceed the local mean density. Caustics in the 
inner halo have very low contrast, however. 



namical mixing, because it is closely related to the number of or- 
bits a particle has executed over its entire dynamical history. We 
examine this using Fig.[T2]which is constructed in exactly the same 
way as Fig. QT] except that the particles in each z = spherical 
shell were ranked by caustic count rather than by stream density. 
Qualitatively the two series of plots show some similarities, but the 
differences between the red and blue populations are substantially 
more marked in the caustic count case. The red points are almost all 
very close to the centres of collapsed clumps at all times, whereas 
the blue population is always diffuse except at low redshift when 
some of it is accreted onto the main halo. The evolution of the blue 
particles in Fig. [12] shows an interesting shell-like structure which 
can be understood with reference to a spherical infall model. Many 
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Figure 14. Local ratio of the spherically averaged intra-stream annihilation 
rate to the spherically averaged smoothed annihilation rate as a function 
of radius. The smoothed annihilation rate is calculated from a SPH-based 
density estimate using 64 neighbours. The intra-stream annihilation rate is 
dominated by contributions around caustics. The standard "boost factor" 
due to unbound small-scale structure is one plus the quantity plotted on the 
vertical axis. Clearly, caustics play essentially no role in enhancing the dark 
matter annihilation luminosity of the inner halo. Around the virial radius 
they contribute about 10% to the annihilation signal. 

of the blue particles have yet to pass a caustic at z = 0, and so are 
on their first passage through the halo. Such particles come from a 
narrow range of Lagrangian radii in the initial conditions, hence the 
apparent shell. The fact that the caustic count is a much better indi- 
cator of the overall level of dynamical mixing than the stream den- 
sity is easily understood; the stream density varies strongly along 
the trajectory of each individual particle reaching high values each 
time it passes a caustic, whereas the caustic count increases mono- 
tonically with the number of orbits completed. 

3.4 Caustics and dark matter annihilation 

The number of fine-grained streams passing through the Solar Sys- 
tem is of interest for the search for dark matter using laboratory 
devices. Caustics, on the other hand, could be of interest for indi- 
rect dark matter searches that focus on the annihilation products 
of dark matter particles. Annihilation explanations for many of the 
apparent "anomalies" in indirect detection signals, for example the 
rise in the positron fraction detected by PAMELA above the ex- 
pectation from secondary production, require annihilation rates to 
be substantially "boosted" above predictions based on "standard" 
particle physics assumptions and haloes with a locally smooth dark 
matter density field. The annihilation rate per unit volume is pro- 
portional to the square of the local dark matter density and so can 
reach very high values in caustics. Caustics could thus, at least in 
principle, provide the necessary boost. To quantify this with our 
simulations, it is important not only to identify caustics, but also to 
estimate their peak densities, since emission around the pe ak dom- 
inates caustic luminosity (see lWhite & Vogelsbergerll2009h . 

In Fig. |T3] we plot median peak caustic density as a function 
of radius for Aq-A. To make this plot we recorded the peak caus- 
tic density and the radial position at which it was achieved for all 
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particles that pass through a caustic in a short time interval around 
z = 0. We bin these caustic passages into a set of radial shells, and 
then find the median value of peak caustic density for each shell. 
The top panel shows results at three different resolutions, Aq-3, Aq- 
A-4 and Aq-A-5. Again, it is clear that the lower resolution of Aq- 
A-5 has significantly affected the results, while Aq-A-3 and Aq-A- 
4 agree reasonably well. The peak caustic density calculation, like 
that of the stream density, depends on the full phase-space distor- 
tion tensor, so it is not surprising that the convergence properties of 
the two measures are similar. In the bottom panel we compare the 
spherically averaged mass density profile of Aq-A-3 with the radial 
distribution of peak caustic density. The thick red line repeats the 
median profile from the top panel, while thin blue lines show the 
upper and lower quartiles of peak caustic density as a function of 
radius. Notice that at each radius the distribution of peak caustic 
densities is very broad, mirroring the range of fine-grained stream 
densities seen in Fig. [7] The peak densities of almost all caustics 
in the inner part of the halo are substantially below the local mean 
halo density. A quarter of all caustics have peak densities exceeding 
the local mean at about 0.4r2oo, but only outside 0.8r2oo do more 
than half of the caustics have peak densities higher than the local 
mean density. 

The low peak densities predicted for caustics in the inner 
halo, suggest that any annihilation boost will be small. This was 
al ready the case for t he isola ted, smoothly growing halo studied 
in lVogelsberger et al.l d2009bl) . In a ACDM halo, additional small- 
scale structure should decrease stream densities, and thus peak 
caustic densities, even further. To investigate this, we calculate the 
ratio of intra-stream annihilation (both particles belonging to the 
same fine-grained stream) to inter-stream annihilation (the two par- 
ticles belonging to different fine-grained streams). The former can 
be calculated for each simulation particle from its GDE-estimated 
stream density and c an be integrated correct l y thro ugh caustics us- 
ing the formalism of IWhite & VogelsbergeJ 1 20091) . The latter can 
be calculated from a SPH kernel estimate of the local smooth dark 
matter density. For the boost to be substantial the ratio of these 
two annihilation rates should be large. In Fig. [14] we plot this ra- 
tio for Aq-A as a function of radius, after averaging the rates over 
thin spherical shells. Again, results for Aq-A-3 and Aq-A-4 agree 
well, while those for Aq-A-5 are clearly affected by its low resolu- 
tion. This reflects the similar effects seen in earlier plots of stream 
and caustic densities. Fig. [14] shows that the contribution of caus- 
tics to the overall annihilation luminosity is very small - boost- 
ing is completely negligible in the inner halo, and is still small 
(about a factor of 1.1) near the virial radius. This is even lower 
than the already sma ll effects seen in the smooth halo collapse sim- 
ulation presented bv lVogelsberger et al.l J2009bl) . Since an annihi- 
lation interpretation of the PAMELA results requires a boost factor 
of 100 to 1000, caustics are clearly far too weak to provide an ex- 
planation. Since bound subhaloes are also unable to produce local 
boosts much above unity near the Sun JSpringel et alj2008bl) . non- 
standard particle physics such as Sommerfeld enhancement (e.g . 
Sommerfeldl Il93ll; l-Hisano et al.l |2004 120051: Icirelli et al.l 12001 
Arkani-Hamed et alJl2009l ; lLattanzi & Silkll2009l) must be invoked 
to explain the PAMELA data through annihilation. 



3.5 Object-to-object scatter 

In the last few subsections we discussed the fine-grained structure 
of Aq-A in considerable detail and studied how it is affected by 
numerical resolution. It is not, of course, clear how representative 
these results are, since they are based on a single Milky Way-like 
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Figure 15. Top panel: Radial profiles of median caustic count similar to 
that in Fig[5]but for all six Aquarius haloes resimulated at resolution level 4 
with a softening of 3.4 kpc. We plot only the region interior to r2oo in order 
to exclude the substantial "noise" at larger radii due to individual massive 
subhaloes. Bottom panel: Ratio of the individual halo caustic count profiles 
to their mean. 



halo. To quantify the scatter expected as a result of variations in 
formation history, environment, etc., it is important to check some 
of these properties for other haloes. The Aquarius Project simu- 
lated six Milky Way-mass haloes at very high resolution, and so 
provides an opportunity to address this issue. The question of how 
representative these particular haloe s are of the full population o f 
similar mass objects is discussed bv lBovlan-Kolchin et al.l d2009l) . 
Here, we concentrate on the caustic count profile which we showed 
in Fig[5]to be very well converged in Aq-A-3, Aq-A-4 and Aq-A-5, 
at least within r2oo where the "noise" due to substructure is small. 
Given this robust result, we decided for the following comparison 
to rerun the other five Aquarius haloes at resolution level 4 in order 
to save computation time. All resimulations were carried out with 
our standard softening length of 3.4 kpc. 

In the top panel of Fig. [15] we show profiles of median caus- 
tic count as a function of radius for all six Aquarius haloes. We 
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plot results within 7-200 only, since at larger radii the profiles are 
subject to large stochastic fluctuations at the positions of massive 
subhaloes (see Fig|5j- The variations between the six haloes are rel- 
atively small and are systematic with radius as can be seen in the 
lower panel of Fig. Q3] where we plot the ratio of each of the indi- 
vidual profiles to their mean. Differences are largest in the central 
regions but still scatter by less than 20% around the mean. There 
is some correlation of the caustic count profile with the radial mass 
density profile. Halo es A and C have the most co ncentrated mass 
density profiles (see iBovlan-Kolchin et alj d2009l) ) and also have 
the highest caustic counts in the inner regions, as might naively be 
expected. Halo E, on the other hand, is the least concentrated of all 
the haloes and yet also has a relatively high caustic count near the 
centre. Clearly, the details of halo assembly history do affect the 
median number of caustic passages significantly. This is not deter- 
mined purely by the local dynamical time of the final halo. 

In Fig.[T6]we do the same exercise for the number of streams. 
For each of our resimulated Aquarius haloes we use the median 
stream density of the particles at each radius to estimate a char- 
acteristic stream number, as in the dashed curves of Fig. [8] which 
demonstrate that good convergence is achieved at resolution level 4. 
We plot the profiles out to 2r2oo in this case; at larger radii stochas- 
tic effects due to massive subhaloes again dominate the variations. 
The top panel shows the profiles themselves, while the lower panel 
plots the ratio of each individual profile to their (geometric) mean. 
Halo-to-halo differences here are larger than for the caustic counts, 
but still relatively modest given the large dynamic range of the ra- 
dial variation. It is interesting that the ranking of haloes here is sim- 
ilar but not identical to that in the caustic count profiles of Fig. 1151 
the most concentrated haloes A and C not only have the largest 
caustic counts in their inner regions, but also the smallest numbers 
of streams. This is unexpected since a naive argument might have 
led one to expect that a larger number of caustics would correspond 
to a larger number of streams. 



4 CONCLUDING REMARKS 

We have demonstrated how N-body simulation techniques can be 
extended to follow fine-grained structure and the associated caus- 
tics during the formation of dark matter haloes from fully general 
ACDM initial conditions. We have shown that, for the large par- 
ticle numbers of our standard experiments, our integrations of the 
geodesic deviation equation (GDE) produce distributions of fine- 
grained stream density and caustic structure which are indepen- 
dent of particle number, hence insensitive to discreteness effects 
such as two-body relaxation. This requires somewhat larger gravi- 
tational softenings than are used in traditional N-body simulations, 
and we show in the Appendix that while caustic count distributions 
are insensitive to the assumed softening, the same is not true for 
the stream density distributions. We have resimulated the six Milky 
Way-mass haloes of the Aquarius Project, allowing us to analyse 
the scatter in fine-grained properties among haloes of similar mass. 

By identifying caustic passages along its trajectory, we are 
able to assign a cumulative caustic count to each simulation par- 
ticle which is robust against numerical noise and serves to indicate 
the extent of dynamical mixing in its phase-space neighborhood. 
This caustic count provides an excellent means to highlight sub- 
haloes and tidal streams in r-v r phase-space plots, because parti- 
cles which became part of condensed structures at early times com- 
plete more orbits, and so pass more caustics than particles which 
remain diffuse until accreted onto the main halo. Filtering the parti- 
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Figure 16. Top panel: Stream number, denned as for the dashed curves of 
Fig. [8] for the six Aquarius haloes at level 4 resolution and with our fiducial 
softening length of 3.4 kpc. We plot out to 2r'2oo since at larger radii these 
profiles are affected by individual massive subhaloes. Bottom panel: Ratio 
of the individual halo stream number profiles to their (geometric) mean. 



cle distribution by caustic count decomposes it into interpenetrating 
components which have experienced different levels of dynamical 
mixing. Particles that have passed no caustic form an almost uni- 
form subcomponent, while particles with large caustic count (e.g. 
> 50) are found only in the dense inner regions of haloes. Particles 
with intermediate count outline the skeleton of the cosmic web. 

Direct dark matter detection experiments are, in principle, sen- 
sitive to the fine-grained structure of the Milky Way's halo at the 
position of the Sun. It is especially important to know if a signifi- 
cant fraction of the local dark matter density could be contributed 
a one or a few fine-grained streams, since each of these would be 
made of particles of a single velocity with negligible dispersion. 
If many streams contribute to the local density and none is domi- 
nant, then a smooth velocity distribution can be assumed. Our sim- 
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ulations show the unexpected result that the distribution of fine- 
grained stream density is almost independent of radius within the 
virialised region of dark matter haloes. Only the low-density tail 
of the distribution appears to extend downwards at small radii. Be- 
cause of the extreme breadth of the distribution, a very large num- 
ber of streams (~ 10 14 ) is predicted at the Solar position, but about 
half of the total local dark matter density is contributed by the 10 6 
most massive streams. The most massive individual stream is ex- 
pected to contribute about 0.1% of the local dark matter density. 
This is potentially important for axion detection experiments which 
have extremely high energy resolution and may be able to detect 
0.1% of the total axion energy in a single "spectral line". Apart, 
from this possibility, dark matter detection experiments may safely 
assume the velocity distribution of the dark matter particles to be 
smooth. 

It has been suggested that the very high local densities associ- 
ated with caustics might significantly enhance dark matter annihi- 
lation rates in dark matter haloes, and thus be of considerable sig- 
nificance for indirect detection experiments. Indeed, for completely 
cold dark matter it can be shown that the contributi on from caus - 
tics is logarithmically divergent and hence dominant dHoganl200ll) . 
This divergence is tamed by the small but finite initial velocity dis- 
persions expected for realistic CDM candidates, and calculations 
for idealised, spherical self-similar mod els indicate quite mod- 
est e nhancements for standard WIMPs jMohavaee & Shandarinl 
l2006h . Our methods allow us to identify all caustics during halo 
formation from fully general ACDM initial conditions. We find that 
they are predicted to make a substantially smaller contribution to 
the local annihilation rate than in the spherical model. This is be- 
cause the more complex orbital structure of realistic haloes results 
in lower densities for typical fine-grained streams and thus to lower 
caustic densities when these streams are folded. The enhancement 
due to caustics near the Sun is predicted to be well below 0.1% and 
so to be completely negligible. Only in the outermost halo does the 
enhancement reach 10%. The standard N-body technique of esti- 
mating annihilati on rates from SPH estim ates of local DM densities 
(see, for example, ISpringel et alj2008bh should therefore be realis- 
tic, and such estimates will not be significantly enhanced by caus- 
tics. Similarly, caustics cannot be invoked to provide the large boost 
factors required by annihilation interpretations of 'anomalies' in re- 
cently measured cosmic ray spectra from the PAMELA or ATIC 
experiments. Given that unresolved smal l subhaloes also appea r in- 
sufficient to provide a substantial boost dSpringel et al.ll2008bh . an 
annihilation interpretation of these signals is only tenable with non- 
standard annihilation cross-sections. 

Our comparison of results for the six different haloes of the 
Aquarius Project showed that the halo-to-halo variation in the fine- 
grained properties we have studied is relatively small, and thus does 
not have significant impact on the applicability of our principal con- 
clusions to the particular case of our own Milky Way. Thus our final 
conclusion must be that the fine-grained structure of ACDM haloes 
has almost no influence on the likelihood of success of direct or 
indirect dark matter detection experiments. 

Although we demonstrated that two-body relaxation and other 
discreteness effects have no significant influence on our conclu- 
sions, we found that our predicted stream density distributions are 
strongly affected by gravitational softening, which modifies the 
tidal forces on particles passing within a softening length or two of 
halo centre. This numerical effect remains a significant uncertainty 
in our results. The analysis in the Appendix shows that reducing 
the softening below our standard value tends to reduce the density 
of streams, thus to suppress further the importance of fine-grained 



structure for detection experiments. However, the change in orbit 
structure induced by our standard softening is similar in scale and 
amplitude to that caused by the neglected baryonic components of 
the galaxies, and even the sign of the effect of the latter is uncertain. 
This is clearly an area where further work is required to reduce the 
remaining uncertainty in our quantitative conclusions. 
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Figure Al. Top panel: Distribution of caustic passages in Aq-A-3 as a func- 
tion of radius and of gravitational softening. Thick lines show the median 
and thin lines the upper and lower 25%, 5% and 1% quantiles of the caustic 
count distribution at each radius. Red lines refer to our standard simulation 
with our fiducial softening length of e = 3.4 kpc. Green and blue lines 
show results for softening lengths half and twice this value, respectively. 
Clearly softening does not significantly affect these distributions. Bottom 
panel: The distribution of stream density in Aq-A-3 (in units of the cosmic 
mean) as a function of radius and of gravitational softening. Continuous 
curves give the median, and dashed and dotted lines the 0.5, 2.5, 10, 90, 
97.5 and 99.5% points of the distribution at each radius. Red curves cor- 
respond to our fiducial softening length of e = 3.4 kpc, green and blue 
lines to softenings half and twice as big, respectively. The stream density 
distribution is significantly affected by gravitational softening at all radii, 
although the effect gets weaker in the outer halo. 



APPENDIX A: EFFECTS OF GRAVITATIONAL 
SOFTENING 

According to Figs.[5]and|7]the caustic count and stream density dis- 
tributions are independent of particle number at our two highest res- 
olutions. This demonstrates that discreteness effects, in particular 
two-body relaxation, are not significantly affecting the phase-space 
quantities we derive from our GDE integrations. In this Appendix 
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Figure A2. Distribution of halo particles at z = in the r/r200 versus 
^min(0.2 — 1) plane for the Aq-A-4 simulation with half our fiducial soft- 
ening, i.e. e = 1.7 kpc. Even at large radius, a surprisingly large fraction 
of particles have minimal pericentric distance within one or two softening 
lengths. 



we investigate whether these distributions are sensitive to our as- 
sumed gravitational force softening. To facilitate this, we keep the 
particle number fixed and change the softening length. Specifically, 
we resimulated Aq-A-3 using softening lengths half and and twice 
our fiducial value, e — 3.4 kpc, but with all other parameters held 
to their original values. 

The top panel of Fig. lATI shows that the caustic count distribu- 
tion is independent of softening in the main halo, but that its high 
count tail is affected by softening in subhaloes. Subhaloes are ef- 
fectively small N systems, and their centres are significantly better 
resolved with a smaller softening length. This results in larger caus- 
tic counts at the centre of subhaloes, which are evident as enhanced 
"noise" in the high-count tail at larger radii for the simulation with 
the smallest softening. This effect is almost absent in the main halo. 
We conclude that our main results for caustic counts are robust to 
ch anges in softening l ength , a conclusion supported by the results 
of IVogelsberger et all J2008h . 

The bottom panel of Fig. lAll shows a similar plot but now for 
the stream density distribution. Here we see a very different result. 
The stream density distribution depends quite strongly on soften- 
ing, with smaller softenings resulting in a shift of the entire stream 
density distribution towards lower values. The effect is largest at 
small radii but is evident throughout the main halo. In the inner 
halo the median stream density drops by more than three orders of 
magnitude when the softening length is reduced by a factor of 4. 
The difference is even larger if we focus on the low-density tail of 
the distribution and is still quite substantial in the high-density tail. 

To understand the origin of this effect, we repeated the Aq-A- 
4 simulation with half the fiducial softening and stored the particle 
data at all ~ 7000 time-steps between z = 4 and z = 0. This al- 
lows us to follow the orbit and the phase-space density evolution 
of each particle in detail. We find the 200 most bound particles at 
z = and use their centre-of mass position to define the halo cen- 
tre at all earlier times. (We checked that they are all indeed still 
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Figure A3. Two-dimensional histograms of properties of a particle sub- 
set defined by the following criteria: halocentric distance at z = be- 
tween 0.01 r2oo and r 2uu ; number of caustic passages from a = 0.26 
to a = 1 between 38 and 42; already part of the main subhalo (i.e. not 
part of any identified subhalo) at z = 2.9. Top panel: Distribution in 
the p Sl min(0.26)/p a , m i n (l) versus r m i n (0.26 — 1) plane. Particles with 
smaller minimum radius over the z = 2.9 to z = period typically de- 
crease their stream density by a larger factor, i.e. Ps.mm(0.26)/p s , m i n (l) 
tends to be larger for smaller r m i n (0.26 — 1). Bottom panel: Distribution in 
the p s m in(0-26)/p a m i n (l) versus p s m i n (0.26) plane. There is no cor- 
relation between the two quantities. Within this particle set, the change in 
stream density after z = 2.9 is independent of a particle's prior history. No- 
tice also that the distribution of the change in stream density is broader than 
the initial distribution, showing that the final stream density distribution is 
determined primarily by effects during the redshift interval considered. The 
small boxes labelled A to D select particles at the edges of the distribution. 
In Fig. lA4l we plot the evolution of stream density and halocentric distance 
for six particles chosen at random from each box. 
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very close to halo centre at 2 = 4.) For every halo particle we 
then construct trajectories of halocentric distance and fine-grained 
stream density, and define a minimum physical halocentric distance 
over any chosen scale factor interval (r m i n (ai — 02)), and a min- 
imum physical stream density achieved between 2 = 4 and any 
later epoch (p s , m i n (a)). Fig. lA2l shows the 2 = distribution of all 
particles in the r/r2oo versus r m i n (0.2 — 1) plane. This histogram 
shows the expected trend that particles at small 2 = radii tend to 
have passed closer to halo centre in their past than particles at large 
2 = radii. Interestingly, the typical minimum radius is more than 
an order of magnitude smaller than the final radius, and the scatter 
in the distribution is quite large. Thus at a third the virial radius 
(around 80 kpc) the distribution peaks at a closest approach dis- 
tance of about 3 kpc, thus within our fiducial softening length. At 
a tenth the virial radius, most particles have passed within 1 .0 kpc 
of halo centre. Thus, everywhere within the halo a significant frac- 
tion of particles have passed through a region where tidal forces are 
substantially affected by the softenings we are using. The tidal field 
sources evolution in the GDE, so this affects our stream density es- 
timates. Specifically, particles feel stronger shearing at pericentre if 
the softening length is smaller, so we expect larger stream density 
decrements for smaller softenings. 

To investigate this further, we select a subset of particles ac- 
cording to the following criteria: halocentric distance at 2 = be- 
tween 0.01 r2oo and r'2oo; caustic count increment from a = 0.26 
to a = 1 between 38 and 42; part of the main halo (i.e. not part of 
any resolved subhalo) at a = 0.26. This selects a set of particles 
with similar orbital periods which have all been part of the main 
body of the system since 2 = 2.9. The top panel of Fig. [A3] shows 
the distribution of these particles in the p s ,mm(0.26)/p s ,mm(l) 
versus r m i n (0.26 — 1) plane. The ratio Ps,mm(0.26)/p s ,mm(l) in- 
dicates the decrease in minimal physical stream density between 
2 = 2.9 and 2 = 0, which, as expected, is strongly correlated 
with minimum halocentric distance. Particles with stream density 
changes substantially smaller than the median have typically never 
been closer to the centre than 2 or 3 kpc, while particles with sub- 
stantially larger than average changes in stream density have al- 
most all been within 3 kpc of halo centre, many within 1 kpc. 
The bottom panel of Fig. I A3 1 shows the same particle set in the 
p s ,min(0.26)/p s , m in(l) versus p s , min (0.26) plane, i.e. the drop in 
physical stream density from 2 = 2.9 to 2 = versus the initial 
stream density at 2 = 2.9. There is no correlation between these 
quantities, showing that the change in stream density after 2 = 2.9 
is independent of what happened to the particle at earlier times. In 
addition, the distribution of stream density changes is broader than 
the initial distribution, showing that the main features of the final 
distribution were established over the time interval considered. 

Together these plots show that, for this particle subset, low 
stream densities at 2 = are associated almost exclusively with 
particles that pass within our fiducial softening length of the centre 
some time after 2 = 2.9. Even in the main part of the stream density 
distribution about half of the particles pass within 3 kpc of halo 
centre. As we will see below, the time-average halocentric distance 
of these particles is about 40 kpc. This explains why changing our 
softening between 1.7 and 6.8 kpc has effects out to large radii, as 
seen in the stream density distributions of Fig. IA1I and why the 
effects are stronger near halo centre and in the low-density tail of 
the distribution. 

The stored data for our test simulation allow us to trace the 
detailed evolution of stream density for individual particles. We il- 
lustrate this for a few particles in the extreme regions of the dis- 
tribution labelled A to D in the bottom panel of Fig. I A3 1 We pick 



6 particles at random from each region and plot in Fig. IA4I their 
physical stream density and physical halocentric distance as a func- 
tion of scale factor. These plots demonstrate that particles in fully 
general CDM halos show quite similar stream density behaviour to 
those in simpler halo models: caustic "spikes" occur several times 
per orbit at turning points of the fundamental oscillations, and sec- 
ular evolution is evident in the steady decrease in the minimum 
stream densities achieved bet ween caustic passages. Th is resem- 
bles the isolated halo results in lVogelsberger et"af d2008h . showing 
that our GDE scheme works correctly in the cosmological setup of 
this paper. The selection criteria for our particle subset are evident 
in the similar apocentric distances and similar caustic and pericen- 
tric passage numbers of all orbits, despite their very different initial 
and final stream densities. 

The stream density decrease is different in the four re- 
gions, reflecting their location in the p s ,mm(0.26)/p s , m in(l) ver- 
sus Ps,min(0.26) plane. Our GDE approach finds stream densities 
spanning almost 18 orders of magnitude for orbits of very similar 
size and period. The orbits in regions A and C look similar both 
in their stream density and in their radial orbit behaviour. The only 
difference is an offset of about 6 orders of magnitude in stream 
density which is preserved throughout the evolution. This empha- 
sises that evolution from 2 = 2.9 to 2 = is independent of what 
happened before 2 = 2.9 if one restricts oneself to main halo par- 
ticles with similar orbital periods. The particles in boxes B and D 
also have similar orbital periods, but they now have similar initial 
stream densities and very different stream density evolution. Those 
in box D barely change their stream density over the period shown, 
while the stream densities of those in box B drop by about 1 1 orders 
of magnitude. Here there are clear qualitative differences in orbital 
structure. The particles in box D avoid the halo centre and in some 
cases are almost circular, and their caustics are relatively uniformly 
spaced in time and show rather little variation. In contrast, the par- 
ticles in box B make much larger radial excursions, often passing 
close to the centre, and their caustic structure is more variable with 
caustic passages often occurring in groups. This suggest that fur- 
ther work on the structure of the stream density evolution curves 
might allow classification of the orbits into tubes, boxes and other 
families. 
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